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ABSTRACT 

The hadronic model of Active Galactic Nuclei and other compact high energy astro- 
physical sources assumes that ultra-relativistic protons, electron-positron pairs and 
photons interact via various hadronic and electromagnetic processes inside a magne- 
tized volume, producing the multiwavelength spectra observed from these sources. A 
less studied property of such systems is that they can exhibit a variety of temporal 
behaviours due to the operation of different feedback mechanisms. We investigate the 
effects of one possible feedback loop, where 7-rays produced by photopion processes 
are being quenched whenever their compactness increases above a critical level. This 
causes a spontaneous creation of soft photons in the system that result in further 
proton cooling and more production of 7-rays, thus making the loop operate. We 
perform an analytical study of a simplified set of equations describing the system, in 
order to investigate the connection of its temporal behaviour with key physical pa- 
rameters. We also perform numerical integration of the full set of kinetic equations 
verifying not only our analytical results but also those of previous numerical studies. 
We find that once the system becomes 'supercritical', it can exhibit either a periodic 
behaviour or a damped oscillatory one leading to a steady state. We briefly point out 
possible implications of such a supercriticality on the parameter values used in Active 
Galactic Nuclei spectral modelling, through an indicative fitting of the VHE emission 
of blazar 3C 279. 



Key words: 

rays: galaxies 



astroparticle physics 
- galaxies: active 



radiation mechanisms: non-thermal - gamma 



1 INTRODUCTION 

The idea that high energy protons can be produced in 
Active Galactic Nuclei (AGN) has been suggested by 
iKazanas fe Ellison] j 19861 ) who considered proton accelera- 
tion by shock waves in th e inner region s of th ese objects. 
Following on this idea, ISikora et al.l l| 19871 ) suggested 
that, in this case, relativistic protons will lose energy in 
photohadronic interactions with the abundant soft photons 
rather than via inelastic collisions with the ambient cold 
proto ns. These i deas were later applied t o the blazar 
jets (|Mannheiml Il993l : iMiicke fc Protherod 120011 ) giving 
rise to what is known today as the hadronic model for 
blazar high energy em i ssion - for a review see iBottcheil 
(2007) and iBottcherl (|2010l ). Similar ideas have been 
appli ed also to Gamma Ray Bursts ([Bottcher fc Dermerl 
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120051 ; iRomero. Christiansen, fc Orellanal 120051 ) in order to 
explain the high energy emission from these objects. The 
hadronic models therefore form a viable alternative to the 
commonly used leptonic ones. 

High energy protons will radiate by synchrotron radi- 
ation, as well as by photopair and photopion while inter- 
acting on any soft photon present. These interactions will 
produce secondaries: electron/positron pairs which are pro- 
duced directly from photopair and via charged pion decay 
from photopion interactions as well as gamma-rays via neu- 
tral pion decay again from photopion; there will be also 
neutrino and neutron production coming as byproducts of 



1 There is a difference between hadronic models describing such 
systems and the corresponding ones used in modelling AGN or 
GRB high energy emission. The physical conditions in compact 
binary systems favour inelastic pp-collisions instead of photo- 
hadronic interactions, which for this reason are neglected. 
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photopion. While neutrinos will escape the source without 
any interactions and so their spectrum will be that at pro- 
duction, the created electron/positron pairs will lose energy 
through synchrotron radiation and inverse Compton scat- 
tering while gamma-rays will be absorbed in photon-photon 
collisions. Therefore in order to calculate the emerging pho- 
ton spectrum one has to follow the evolution of these sec- 
ondaries which can be complicated due to the formation 
of intense electromagnetic (EM) cascades initiated, e.g., by 
7— rays from ty° decay. 

Usually the modelling of hadronic processes assumes 
that the target photons come either from an external source 
or from synchrotron radiation of a co-accelerated leptonic 
component. One largely overlooked aspect is the possibil- 
ity that protons interact with their own radiation, for ex- 
ample, with the soft photons produced from the aforemen- 
tioned EM cascades. First a ttempts to incorporate these 
into the models were made by Stern fc Svenssor] l|l99ll ) and 
IStern. Sikora. fc Svensson] (| 1992T ) . Using Monte Carlo sim- 
ulations the above authors found that the system of pro- 
tons and photons can exhibit limit cycles. However, this os- 
cillating behaviour of the system could not be interpreted 
as the result of a specific feedback mechanism, let alone 
studied in a systematic way. One of th e operating feedback 
proces ses was studied analytically by Kirk fc Mastichiadis! 
i|l992t ) using the kinetic equation approacfQ. They showed 
that a sufficient number density of protons can make the 
system unstable, causing runaway pair production: syn- 
chrotron photons of the relativistic electron-positron pairs 
become targets for the protons which produce more pairs. 
The feedback leads eventually to fast proton energy losses. 
This amount of energy lost by the protons in a small 
time interval is transfered to photons and it is seen as 
a flaring event. The operation of this feedback loop was 
late r confirmed numerically by Ma stichiad is fc Kirkl (| 19951 ) 



id lMastichiadis. Protheroe. fc Kirkl ( 20051 ) . who also found 
cases where the system showed a limit cycle behaviour. 

In the present paper we examine another type of feed- 
back which can operate in hadronic systems. For this we 
capitalize on the ide a s of non-linear photon quenchin g 
l|Stawarz fc Kirkl 120071 ; iPetropoulou fc Mastichiadis! l201ll ). 
According to this, 7— rays produced in a spherical volume 
cannot exceed a critical luminosity that depends only on 
the source's magnetic field and radius. If they do, then soft 
photons will be produced automatically and quench the 'ex- 
cessive' 7— rays. In the context of a hadronic system 7— rays 
can either be produced directly by proton synchrotron radi- 
ation or indirectly by photopair and photopion interactions. 
Then the following loop suggests itself: 

(i) Protons cool on soft photons producing 7— rays. 

(ii) The 7— ray luminosity is quenched and turned spon- 
taneously into soft photons which feedback on 1. 

Clearly this mechanism can tap energy stored in protons 
and transfer it to radiation. At the same time it shows that 
the hadronic system is a dynamical one and its behaviour 
can be more complex than it is customarily assumed. A de- 



tailed study of the aforementioned feedback mechanism, in 
the case where 7-rays are the byproduct of photopion inter- 
actions, will be the subject of the present work. The paper 
is structured as follows: In §2 we describe qualitatively the 
system and define two regimes of operation. Next we con- 
struct a system of non-linear equations which we solve first 
analytically in a simplified form (§3). In §4 we show in a 
semi-analytical way the role of various processes on the dy- 
namical behaviour of the system, while in §5 we back our 
results presenting a full numerical study of the problem. 
In §6 we present an indicative astrophysical example, where 
some of our results are applied to the blazar 3C 279. Finally, 
we conclude in §7 with a summary and discussion. 



2 QUALITATIVE DESCRIPTION OF THE 
PHYSICAL SYSTEM 

2.1 Linear regime 

We assume a spherical source of radius R with embedded 
magnetic field B and a monoenergetic proton distribution 
of number density n p and Lorentz factor 7 P . This region is 
also filled with monoenergetic radiation of energy e (nor- 
malized to electron rest mass energy) and number density 
« C >EI, which we will assume comes from outside of the source. 
Thus, we denote it as 'external'. 

The injected high energy protons will interact with the 
external photons through inelastic photopair and photopion 
collisions, provided that the respective threshold conditions 
are satisfied. We will assume that the condition 



£o7p 



(1) 



where m n is the pion mass, is always satisfied. This means 
that both photopair and photopion operate; however since 
the target photons are monoenergetic, it guarantees that 
photopion will be the main loss mechanism for p rotons 
jSikoraet al.lll987l : iBegelman. Rudak. fc Sikoralll990l ). The 
produced charged and neutral pions will decay producing 
electron/positron pairs (for brevity we will refer to them 
simply as 'electrons') and 7— rays. The former will radi- 
ate photons mainly through the synchrotron process, since 
inverse Compton scattering will be greatly suppressed by 
Klein-Nishina effects. We will refer to these synchrotron pho- 
tons as 'hard', since for proton energies with values typical 
of the AGN hadronic models, these can in principle emerge 
in the 7— ray regime. 

One can quantify the above by noting that the sec- 
ondary electrons from the charged pion decay are produced 
with a Lorentz factor 



(2) 



where rj p — k p /4 ~ 0.08 and k p is the inelasticity of the 
interaction assumed to be ~ 0.3. The factor 1/4 arises from 
the assumed energy equipartition between the lepton and 
the three neutrinos produced by the charg ed pion decay - 
see lDimitrakoudis et alj |to appear in 2013 ). Assuming that 



2 A discussion abo ut the different nu merical approaches em- 
ployed is presented in lStern et al. (1995). 



3 We note that here and through the present work tilted quanti- 
ties denote quantities with dimensions. 
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these electrons emit at the critical synchrotron energy 



: &7e,7T, 



(3) 

where b = B/B cr and B cr = 4.413 x 10 13 G the critical 
value of the magnetic field strength, we find that for typical 
values of 7 P = fO 8 and B = 1 G, eh is in the TeV regime. 

Let Etot be the energy loss rate of all protons of energy 
7 P due to interactions with the photons. This energy is dis- 
tributed to the produced secondaries. Assuming that their 
cooling is fast - an assumption which is reasonable since 
both B and 7 C are assumed to have high values, one can 
argue that the energy injected into secondary electrons will 
be instantaneously radiated as hard photons. Then we can 
define the injected hard photon compactness as 



-Etot 0T 



(4) 



h ^ 4nRm cC z ' 

where or is the Thomson cross section and is the fraction 
of energy that goes to the secondary electrons. Furthermore, 

we can connect E^Qt to the single proton energy loss rate 

E p through the relation 



Etot = hpVE p , 



(•>) 



where V is the volume of the source. Since photohadronic 
losses can be considered catastrophic, i.e., a relativistic pro- 
ton can lose a substantial amount of its energy in one colli- 
sion with a photon, we can write 



E p ~ kpEpC / Ax 5p 1 (-ypx)np h {x), 



(6) 



where E p — 7 P m p c 2 is the proton energy, a P7 is the relevant 
cross section, h p h is the target photon population and x the 
target photon energy in units of m c c 2 . Under our assump- 
tions, i.e., catastrophic energy proton losses and monoener- 
getic particle distributions, we can adopt working, from this 
point on, with k p — 1 without loss of generality. Further- 
more, we can approximate the cross section with a Heaviside 
function of the form 



(7) 



with <7p 7 = 10 4 ; for a plot of the total ex pression of the 
cross section see Fig. 3 in iMiicke et al1(|2000h . Using also the 
fact that the only target photons present are the external 
ones, eq.© becomes 



Ep ~ 7 P m p c a pl cn ex . 



(8) 

Combining relations l[4" )l - l[8 jl one can immediately deduce 
that the compactness (or luminosity) of the hard photons 
depends on both n p and h cx . In this case the system can be 
considered to operate in the linear regime, since all cooling 
is provided by the external photons. 



2.2 Non-linear regime 

The previous results indicate that for sufficiently high val- 
ues of hp or n cx , the injected hard photon compactness 
ca n take high values as well. However, as it was shown 
in IStawarz fc Kirk! (|2007l ) and IPetropoulou fc Mastichiadisl 
|201ll ) - henceforth SK07 and PM11 respectively, if the hard 
photon compactness is larger than some critical value 



that depends only on eh and on source parameters such 
as B and R, even small initial perturbations of low energy 
photons present in the source, can grow and lead to an au- 
tomatic quenching of the hard photons. This is a purely 
non-linear process. In this case, electron-positron pairs grow 
spontaneously in the source and the 'excessive' hard radia- 
tion is absorbed by the synchrotron photons emitted by the 
pairs. Thus, a soft photon population of number density n s 
and energy e s appears spontaneously in the source. Assum- 
ing equipartition of energy between the created pairs one 
finds that e s is given by 



&7e 



(SO 



These automatically produced photons have the same en- 
ergy with those produced by the absorption of hard pho- 
tons on the external ones; note that this is a linear process. 
Thus, both linear and non-linear absorption of hard photons 
result in the formation of a third photon population in the 
system with compactness £ s . This new component will start 
playing a role in proton cooling through eq. ([6]), since now 
hph ~ n ax + h s . If its number density grows sufficiently high, 
then it is possible that the relativistic protons will start cool- 
ing more efficiently on them than on the external photons. 
In this case, proton cooling becomes non-linear. 

Figure[T]summarizes the different processes operating in 
the system. Arrows leading to the circles of Fig. [T] denote in- 
jection of the corresponding particle species into the source, 
whereas arrows coming out of the circle of hard photons im- 
ply their subsequent absorption. In order to emphasize the 
existence of two absorbing channels for the 7-rays, the non- 
linear one is shown with a dashed line. We note also that 
the secondary electrons that are the intermediate products 
of the different processes operating in the system and re- 
sponsible for the emission of hard and soft photons are not 
shown in Fig. [1] 

We would like to examine next under which conditions 
the non-linear loop operates. There are two conditions on the 
energies, that must be simultaneously satisfied. The first is 
the feedback criterion for photon quenching. This can be de- 
rived from the requirement that the magnetic field is strong 
enough so that the synchrotron photons of the produced 
pairs are above the threshold for photon-photon absorption 
on the hard photons. This leads to the condition for the 
magnetic field in the source b > Se^ 3 (SK07; PM11). The 
last relation combined with eqs. @, Q and @ sets a lower 
limit to the magnetic field strength which depends only on 
7p 



a/8 ( mc 



1/2 



7 P 



3/2 



(10) 



The second is that the energy of the soft photons is high 
enough for the production of pions in interactions with the 
protons, i.e., the relation 



e s 7p 



(11) 



should hold. This sets another lower limit to the magnetic 
field strength given by 



b>b v 



. 1/3, 

uitt \ I m 



4/3 

P \ -5/3 

Vp 7 P 



(12) 
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It is interesting to note that both b q and b n depend 
only on 7 P . Clearly, in order for non-linearity to appear in 
the system, the (normalized) magnetic field of the source 
should satisfy the condition 



b is max(6 q , b^) 



(13) 



This is not a strict limit. For instance, if 7 P = 10 one finds 
that B q = 0.04 G and = 0.03 G. Thus, magnetic fields of 



the order of 1 Gauss can easily satisfy condition (|13[) . Per- 
haps more limiting are other various effects, that we proceed 
to discuss next: 

(i) Photon quenching is based on the premise that soft 
photons start building in the system once the hard photons 
are above a certain critical luminosity, . The idea is that 
hard photons are absorbed on the automatically created 
soft photons; the produced pairs emit more soft photons 
through synchrotron radiation which causes more pair pro- 
duction on the hard photons etc. However, in the present 
situation, the existence of an external photon population 
complicates the picture as the hard photons might pair 
produce on them, in parallel to the internally built soft 
photon population. This can act as a stabilizing factor, as 
it can inhibit the soft photon build-up in the system. This 
complication can be avoided if one assumes that the relation 
£h£o < 2 holds, i.e., that collisions between the two photon 
populations are below the threshold for pair-production. 
Since only eh depends on B and 7 P , one can find values for 
one of these parameters, in order ehe < 2 to hold (see also 
§3). In the analytical treatment of the next section we will 
first start with this assumption but later we will relax this 
and study the effects, that the hard photon pair produc- 
tion on the external ones has on the dynamics of the system. 

(ii) As the soft photons start building up, the secondary 
electron-positron pairs that produce them will start losing 
energy gradually through inverse Compton scattering 
instead of synchrotron radiation. This means that the 
photon quenching loop could become less efficient, because 
not all of the secondary electron luminosity will end up in 
the energy bin e s . This effect might become important near 
saturation, i.e., when the soft photon luminosity reaches 
its maximum value. As a first step, we will ignore inverse 
Compton scattering in our analytical treatment. Then, we 
will include it in an approximate manner to determine its 
effects. Finally, we will take this fully into account in our 
numerical treatment in §5. 

(iii) The exact nature of the external photon density 
in the source was of no importance so far. As a first 
order approximation to the problem, we neglect proton 
synchrotron radiation, assuming that the initial photon 
distribution is purely external. Under this assumption, both 
the external photon number density n ex and energy e are 
treated as free parameters, instead of being determined by 
source properties, as the magnetic field B and the proton 
energy 7 P m p c 2 . However, in §5, where numerical solutions 
of the full problem are presented we replace the external 
source of photons by the proton synchrotron radiation. 

(iv) In the analysis above the role of the hard photons was 
given to the synchrotron photons of the charged pion decay. 



Protons 
n 



External 
photons 





Figure 1. Schematic diagram of the operating loop of processes 
between protons and photons. 



Without loss of generality we could also assign them to tv° 
initiated secondaries: Assuming that there is equipartition 
between the produced 7— rays from a it — decay and that the 
inelasticity parameter is about the same as that for charged 
pion decay, we can write in full analogy to eq. Q 



2r? P 7 P - 



(14) 



These are extremely hard 7-rays and will always be above 
the threshold for pair production on the external photons 



e Q . The produced pairs will have energy 7 e 



y /2 that 



is exactly the energy of the injected pairs through charged 
pion decay (c.f. eq. ([2])). Therefore, both charged and neutral 
pions decay and produce pairs of the same energy, that cool 
by synchrotron providing the hard photon emission. 



3 ANALYTICAL APPROACH 
3.1 Simplified equations 

In the most general case, the hadronic system consists of 
three species of particles, namely protons, electrons and pho- 
tons. In order to be descr ibed, the kinetic equations for the 
particles must be solved (Mastich iadis fc Kirkj|l995T ). Their 
generic form is: 



Hi 



— ^— H — — — = L l + 1 



(15) 
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where the index i can be one of the subscripts 'p','e' or 
'7' referring to protons, electrons and photons respectively. 
The operators 12 and Q ! denote the losses and injection 
terms respectively, whereas rn/U^sc is just the escape term 
from the source, with each species having its own escape 
time ti, esc- For photons the relation t ltCSC = t cr — R/c 
holds, whereas for protons we adopt t p , CS c = W 3 t CI through- 
out the present work. T he explicit expression s of th e op- 
erators can be found in iMas tichiadi s fc Kirkl (|l995l ) and 
iMastichiadis et al.l (|2005h - henceforth MK95 and MPK05 
respectively. 

The unknown functions to be determined are the num- 
ber densities rk, that can be normalized as follows: 



n P (E p ,T) 



n e (E e ,r) 



with 7 P 



E a 



hr With 7 = 



,(e 7 ,r) = 



GTRm c c 2 



n 7 (e 7 ,r) 



with e-y — 



(16) 



(17) 



(18) 



(Tt-RtTIcC 2 

Time r is normalized to the crossing time of the source t cr , 
i.e., t — t/t CI . From this point on, we adopt working with 
dimensionless quantities. 

For the purposes of an analytical treatment, we make 
two major simplifications: 

(i) Since electron cooling can be considered fast for 
typical values of the system's parameters (see also §2), we 
can neglect the equation of the electrons. 

(ii) We simplify the equations describing the physical sys- 
tem to such a point as to retain only the key processes. 

Thus, we assume a constant monoenergetic injection of pro- 
tons Q P o(7p) into a spherical source of radius R and em- 
bedded magnetic field B. We assume also a monoenergetic 
distribution of external photons n cx at energy e = 7 p ~ 1 - 2i!L , 
that acts as a target for high energy protons. We note that 
initially, n cx is the only distribution of low energy photons 
present in the source. Charged pions produced by proton- 
photon pion processes, decay into electron-positron pairs 
with Lorentz factor 7 CjW (see eq. ([2}), that emit through syn- 
chrotron radiation, hard photons with corresponding num- 
ber density nh(eh)- In order to ensure that inverse Compton 
scattering of the external photons by the aforementioned 
electrons is not as important as synchrotron cooling, and 
therefore can be safely neglected, we assume that u ox << «b 
or equivalently i ax « £b, where m and ii are the energy 
densities and compactnesses respectively. At this point, it is 
useful to define the different compactnesses that appear in 
the present work: 

~3~' 

7pW p 



ii = 



■■ s, h, ex 



l B = °tR- 



II B 



(19) 
(20) 
(21) 



In our treatment, secondary electrons do not affect the dy- 
namics of the system. They rather play an intermediate role 
for transferring energy from hard photons to lower energy 
photons. We examine separately the hard and soft photon 



populations, rih and n s respectively, by writing a kinetic 
equation for each one of them. Thus, the quantities to be 
determined now, are n p , and n B . The physical processes 
to be included into the simplified version of equations are: 

(i) Constant proton injection Q po and proton escape 
Ease — ~ n p/ r p> that act as a source and a loss term re- 
spectively in the proton equation. 

(ii) Proton-photon pion production, that acts as a loss 
term for protons L^_, pn and as an injection term for hard 
photons Qp 7 ^ p ^. 

(iii) Photon-photon pair production, that acts as an ab- 
sorption term for hard photons and as an injection term 
for soft photons Q 77 . 

(iv) Photon escape from the source in a crossing time, 
i.e., L2 SC = Tip- 

Under these considerations, the simplified kinetic equations 
for each species are given by: 

^ rip 



IP. i r p 



-n h + Q 



-n s + 1 



_1_ r n 
p 7 — ►pTr ~ 7"y 



(22) 



(23) 
(24) 



The analytical study of the dynamics of the above set of 
equations can be simplified even further, if we assume that 
hard photons cannot be absorbed by the external ones (see 
point (i) of §2). To ensure this, we use for a given 7 P , such 
values of B that do not allow further absorption of hard 
photons. Thus, the condition eh£o < 2 should hold. This, 
combined with eq. ((2} and sets an upper limit for the 
magnetic field strength 



b < b a = 2- 



(25) 



Equations (|10p . (|12[) and (|25f) can be combined in order to 
create a parameter space of allowed values of B for different 
proton energies - see Fig. [5] Solid line shows the upper limit 
of eq. (|25[) . whereas dashed line shows the maximum value 
of the two lower limits given by eq. (|13|) . 

According to §2, if the compactness of hard photons is 
below the critical value the instability that leads to an 
'automatic' quenching of hard photons cannot grow. Thus, 
there is no injection of soft photons in the system and hard 
photons do not suffer any absorption. This situation corre- 
sponds to the linear regime described in §2.1. In this case, 
assuming that proton losses are catastrophic we can write 

Epry—>p n — £ 7 p 7 ^p7^cx and Qp 7 — >p7r — ^.^p^cx- 

The set of 

equations (|22|) -(I24 | ) degenerates into a system of two equa- 
tions (SI): 



iih 



= — n h + An p n e 



(26) 



(27) 



The normalization constant A of the injection term in the 
hard photon equation above, is found after assuming that 
a fraction ^ of the total proton energy goes to secondary 
electrons, that radiate further all their energy through syn- 
chrotron producing the hard photons. Under these consid- 



erations, one finds that A = 



The steady state 
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where 




Figure 2. Upper (solid line) and lower (dashed line) limit of the 
magnetic field, where the lower limit is given by max(B q , B^) (see 
eq, lll3l l). Any value of the magnetic field that lies below the solid 
and above the dashed line, does not allow absorption of hard pho- 
tons by the external ones. The lines are drawn for the minimum 
energy of external photons that satisfies the energy threshold for 
the photopion process, i.e. e a = m7T . 



solution of system SI is 



An c 



(28) 



(29) 



Two limiting cases can be implied 



where G p — — + <Tp 7 n c 
by the form of G p : 

. Proton escape is more significant than proton cooling 
on the external photons, i.e., <7p 7 n cx <C 1/t p . In this limit, 
the steady state solution for hard photons is proportional 
to the product Q po n ex . One can find a combination of 
values for Q po and n ex that lead to £h > £^ . For example if 
proton cooling is not efficient in injecting hard photons into 
the system because of a low n ex , a high value of the proton 
injection rate is needed and vice versa. 

. Proton cooling on the external photons is more signif- 
icant than proton escape, i.e., <7p 7 n cx > l/r p . In this limit, 
n| s oc Qp . Thus, the system can become supercritical for a 
high enough value of the proton injection rate. 

So far, the proton-photon system operates in a linear 
'subcritical' regime with a very well described behaviour. 
In order to investigate how the non-linear terms affect the 
evolution of the system, we will focus only on cases where 
quenching is relevant. If hard photons are being injected into 
the system with t^ s > (.^ , then a soft photon population 
n s (e s ) appears because of the quenching of the 7-rays and 
the system becomes 'supercritical'. In this case an additional 
equation for the soft photons is required. Thus, the set of 
equations (|26[) - (|27p becomes (system S2): 



= — nh + An p n cx + An p n E — Chn s rih 
= — n B + C s n s rih, 



(30) 



(31) 
(32) 



and 



C 3 = 



'77 

*2 ' 



(33) 



The term <r 77 that appears in expressions ([33} is given by 



eh£scr 77 (ehes) 



(34) 



with the exponent 's' denoting the absorption on the soft 
photons and <r 77 the cross section of photon-photon absorp- 
tion, measured in units of the Thomson cross section gt- 
Here and t hroughout this work , we us e the approximate ex- 
pression of ICoppi fc Blandfordl l)l99(t l: 



cr 77 (a;) = 0.652 



(* 2 - 1) 



\n(x)H(x - 



(35) 



where x is the product of the dimensionless photon energies 
and H{x) is the Heaviside step function. Since in our analysis 
we have assumed monoenergetic photon distributions and a 
cross section approximated by a step function, a 77 in defi- 
nitions (|33p is just a constant, that takes different values for 
different pairs of photons. The last terms in the right hand 
side of eqs. (|31[) and (|32[) account for the photon-photon 
absorption. The created pairs are the intermediate products 
that transport their energy through synchrotron radiation to 
soft photons. The constants Ch and C s given above are de- 
termined by energy conservation. As far as the soft photons 
produced have enough energy to satisfy the energy threshold 
for photopion interactions, (see relation (HU), they can act 
as an additional target for protons. This explains the last 
term in the right hand side of eq. (130[) . This additional sink 
term of protons has a corresponding injection term An p n s , 
that appears in eq. (|31fl . 

We note that an expression of the critical compactness 
for the automatic quenching of hard photons has been pre- 
sented by SK07 and PM11. In the present analysis we have 
made certain simplifying assumptions that differ from those 
in the aforementioned papers. Thus, for reasons of consis- 
tency, we calculate the modified expression of 



(I 



6h 

3a 



(36) 



(see Appendix C for the derivation). 



3.2 Dynamical study of the system 

After having set the framework of the physical problem 
and clarified our assumptions, we proceed to investigate the 
mathematical properties of the system. 

Setting in equations (I30|l - (|32|l the time derivatives equal 
to zero we can determine the possible steady states, that ac- 
tually are the 'fixed points' of the dynamical system. The 
system of equations S2 has one or three fixed points depend- 
ing on the existence or not of the coupling terms between 
hard and soft photons. We will treat the system of equations 
that obtains three fixed points, since this can in principle 
show non-linear temporal behaviour. 

The first fixed point Pi is just the steady state solution 
presented in the previous section, i.e., Pi (n^ s ,0,np s ). Per- 
turbations of this steady state solution can either grow or 
decay with time as e ET . In the first case, the fixed point Pi 
is unstable whereas in the second case it is stable (see Ap- 
pendix A for a detailed analysis). The growth or decay rate 
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s of the perturbations as a function of the proton injection 
rate, for a fixed external number density, is shown in Fig. [3] 
It is interesting to note that the exponent s for the growing 
solutions depends strongly on Q po , whereas the respective 
one for the decaying solutions is rather insensitive to Q po 
and close to —1, denoting the free escape of the produced 
soft photons in a crossing time. The star denotes the critical 
value of the proton injection rate, Q po , above which s > 
and an automatic build-up of soft photons in the system is 
possible. This value depends on the external photon density 
as: 



Qpo(^c 



C B An c 



(37) 



Low values of rt cx suggest inefficient proton cooling and 
therefore inefficient injection of hard photons - see eqs. (1261) - 
(|27p . Only if the injection rate of protons in the source is 
high, can the hard photon compactness increase sufficiently 
leading the system to instability. In other words, for a given 
rt cx one can find always a sufficiently high Q po in order to 
make the system unstable. Qp r D is directly related to a criti- 
cal proton number density n p given by 



(38) 



1/Yp + <r° 7 n e 



Thus, the problem of choosing a suitable pair of values 
(<5 P o, n ex ) that lead the system to instability, can be trans- 
ferred to the problem of loading the source with a critical 
proton content. Let us now consider the inverse of the above 
statement: if the injection rate of protons is given, is there 
any value ofn ex that will make the system unstable? In this 
case, we find that the external number density must satisfy 
the following condition: 



CsAQ p 



(39) 



with the constraint Q po > -^j. Thus, for low enough pro- 
ton injection rates, the system cannot become unstable. The 
above remarks lead to the conclusion that, between the two 
parameters, Q po and n ex , the fundamental one in making 
the system unstable is the proton content of the source. 

When the initial growth of the perturbations is en- 
sured, the subsequent behaviour of the system depends on 
the properties of the second fixed point P2. The system 
can either reach the steady state P2 or vary periodically, 
making a limit cycle in phase space. This behaviour can 
be predicted by calculating the eigenvalues of the matrix 
Mi of the linearized system of equations, near the point 
P2. For different values of the physical parameters of the 
problem, we find always one real negative eigenvalue, Ai. 
The two remaining eigenvalues A2,3 can either be com- 
plex conjugates or both real and negative (see Appendix 
B for more details). Table 1 summarizes the dynamical 
behaviour of the system. We note that, for the terminol- 
ogy regarding the classification of a fixed point of a three- 
dimensional system, we have adopted th e one presented in 
iTheisel. Weinkauf. Hege. fc Seidell (120031 ). 




Figure 3. Growth/decay rate of the perturbed proton and pho- 
ton densities from the steady state PI as a function of the proton 
injection rate. Above a certain value of Q po , that is denoted with 
a star, only growing perturbations exist. Parameters used for this 



plot are: B = 0.7 G, 7 P = 2 X 10 7 , t a 



: 7 P 



and n c 



Eigenvalues 


Classification of point P2 


Dynamical behaviour 


A, < 0, for i= 1,3 


attracting node 


steady state 


Ai < 0, A 2 = X* 


focus 




Re(A 2 ) = Rc(A 3 ) > 


repelling saddle 


limit cycles 


Re(A 2 ) - Rc(A 3 ) < 


attracting 


damped oscillations 



Table 1: Classification of the fixed point P2 and of the ex- 
pected dynamical behaviour of the system, based on an eigen- 
value/eigenvector analysis of the corresponding matrix M2. 

It is interesting to investigate how the two main pa- 
rameters of the physical problem, i.e., the proton injection 
rate Q po and the external number density n ex , are related 
to the dynamical behaviour of the system. For this purpose, 
we calculate the eigenvalues of matrix M2 around the fixed 
point P2 for different values of Q po and n cx , having first en- 
sured that the fixed point of the steady state Pi is unstable 
(see Appendix A). The results are presented in the next two 
paragraphs. 



3.2.1 Dependence on proton injection rate 

We assume first that the external density n cx is fixed to a 
certain value and study the effects of the injection rate Q po . 
Figure [4] shows the calculated eigenvalues \2,s for different 
values of Q po . The conclusions drawn from it can be sum- 
marized in the following points: 

(i) Starting with low values of Q po , we find that A3 = A2. 
As long as Re(A) > 0, point P2 acts as a repelling focus - 
saddle point. Having also ensured that the fixed point of 
the steady state Pi is unstable, the system follows in phase 
space a closed periodic trajectory, i.e., a limit cycle that is 



© ... RAS, MNRAS 000,[TJfT9] 



8 M. Petropoulou and A. Mastichiadis 



moreover stable. 



(ii) For higher values of the proton injection rate we 
find complex eigenvalues with Re(A) < 0. Therefore, the 
point P2 acts as an attracting focus - saddle point. In 
this case, the system settles down in a new steady state 
(given by the fixed point Pi). In phase space, this cor- 
responds to a 'spiraling' trajectory ending at the fixed point. 

(iii) Finally, for even higher values of Q po both eigenval- 
ues become real and negative. Point P2 can be characterized 
as an attracting node. The physical system reaches the new 
steady state very fast, showing no oscillations. 

An example of such a transition in the dynamics of the sys- 
tem, is shown in figures [5] and [6] The solution shown with 
dashed-dotted line corresponds to a limit cycle case with 
period ~ 170 i C r- For reasons of clarity, Fig. [S] zooms in the 
early time behaviour of the system. Therefore, only the first 
and a half cycle of the periodic solution is shown. The solu- 
tions presented in the aforementioned figures are found after 
integrating the system of equations S2 with initial conditions 
w-h(O) = n p (0) = and n s (0) = e, with e — > indicating an 
initial perturbation of soft photons in the source. The pa- 
rameters used ensure that at some point the compactness of 
hard photons becomes larger than fjf . Each of the values of 
Qpo used in this example corresponds to a different point of 
Fig. [4] The conclusion is that the temporal behaviour of the 
system is very sensitive with respect to the proton injection 
rate. 

For reasons of mathematical consistency we make also 
the following remark: As the proton injection rate increases, 
the real part of the complex eigenvalues changes sign and 
from positive it turns into negative. At this transition a sec- 
ond unstable limit cycl^3 surrounding the fixed point P2 ap- 
pears inside the stable limit cycle (see point (1) above). An 
integration of the equations of S2 with initial conditions ly- 
ing outside the stable limit cycle would show that the system 
falls onto this cycle instead of 'spiraling' down to the fixed 
point, as one would expect according to point (2). As the 
varying parameter Q po increases further, |Re(A2,3)| increases 
too and the unstable limit cycle approaches the stable one. 
At some value, which for the specific example of Fig. [3] is 
log Qpo = —10.9, the two limit cycles coalesce and the sys- 
tem undergoes a bifurcatiorjf] since the phase space changes 
qualitatively. From this moment on the system can fall into 
the fixed point P%, as described in (2) above. Thus, the con- 
dition Re(A2,3) < that occurs for logQ po = —11.5 in our 
example, does not ensure strictly speaking the damped os- 
cillatory behaviour of the system. 

It is beyond the scope of the present work to proceed 
into a detailed study of the bifurcation mentioned earlier, 
since for the initial conditions of physical relevance to our 
analysis, the existence of the bifurcation does not affect the 




log Q 



Figure 4. Plot of the two eigenvalues A2,3 for 
function of the injected proton rate. As long as the eigenvalues 
are complex conjugates, their real and imaginary parts are shown 
with open and filled circles respectively. At a certain value of 
Qpo both eigenvalues become real and negative (open and filled 
diamonds). The solid line corresponds to the null value. The other 
parameters used are the same as in Fig. \3\ 
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Figure 5. Time evolution of the proton distribution for different 
values of the proton injection rate logQ po = —11.15 (dashed- 
dotted line), logQpo = -10.8 (dotted line), logQ po = -10.2 
(dashed line) and logQ po = —9.2 (solid line). The initial con- 
ditions for each numerical run are r»h(0) = ip(0) = and 
n s (0) = e — > 0. The number density of the external photon distri- 
bution is n cx = 2 in all cases. All the other parameters used are 
the same as in Fig. [3] 



qualitative features of the transition from the limit cycle 
phase to the damped oscillatory one. 



4 We note that the existence of the unstable limit cycle was found 
numerically while integrating the equations of system S2 and not 
analytically. 

5 The existence of the bifurcation is not a characteristic property 
of the system for all values of the parameters. For example, if 
n cx = 10 no bifurcation of this type is found. 



3.2.2 Dependence on number density of external photons 

In this section we show an example of the effects of the exter- 
nal number density n ex on the dynamical properties of the 
system. An increase of the external number density leads 
to the same transition of the dynamical behaviour of the 
system, as the one described in the previous section. Figure 
[7] shows the calculated eigenvalues for a large range of n ex 
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Figure 6. Two-dimensional plane of the phase space for different 
values of the proton injection rate logQ po = —11.15 (dashed- 
dotted line), logQ po = —10.8 (dashed line), logQ po = —10.2 
(dotted line) and logQ po = —9.2 (solid line). Same parameters 
used as in Fig. \E\ 



Figure 7. Plot of the two eigenvalues A2,3 as a function of the 
external number density n ex for injected proton rate log Q po = 
— 11.15. While the eigenvalues are complex conjugates, their real 
and imaginary parts are shown with open and filled circles respec- 
tively. At a certain value of n ex both eigenvalues become real and 
negative (open and filled diamonds). The solid line corresponds 
to the null value. All other parameters used are the same as in 
Fig. [3] 



extending up to high values. The qualitative features of this 
plot are the same as those of Fig. [4] One should keep in mind 
though, that the conclusions drawn from Fig. [7J regarding 
the dynamics of the system are not valid for the whole range 
of values of n ex shown. The reason is the following: As n ex 
increases, the compactness of the external photons / cx in- 
creases too. At some value, which for the specific example 
is just n ex = 35, it becomes larger than the compactness of 
the magnetic field Is = 1.6 x 10 -4 . This means that the sec- 
ondary electrons produced by the pion decay, cool preferably 
through inverse Compton scattering on the external photon 
field, rather than through synchrotron radiation. Thus, the 
system of equations S2 we used to make our mathematical 
analysis is physically not valid. However, aim of Fig. [7] is to 
simply show the mathematical similiraties of this case with 
the previous one. 

Figure [8] shows a two-dimensional plane of the phase 
space for two different values of the external density. All 
the other parameters are kept constant. The corresponding 
time evolution of the proton distribution is shown in Fig. [9] 
External photons act as a stabilizing factor for the system, 
since high enough values lead the system to a steady state. 



4 ENHANCING NON-LINEARITY WITH 
ADDITIONAL PROCESSES 

4.1 Photon-photon absorption by the external 
source 

In the previous section, in order to avoid the appearance 
of more terms in the equations and be able to treat them 
analytically, we restricted our analysis to values of B taken 
from the parameter space shown in Fig. [2] and to exter- 
nal photons having energy equal to the threshold energy for 
proton-photon pion processes. If one wishes to use values 
of the magnetic field more related to astrophysical sources 
or consider more energetic external photons, the absorption 
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Figure 8. Two dimensional plane of the phase space for two 
values of the external number density n ex = 2 (solid line) and 
n cx = 4 (dashed line). The proton injection rate is logQ po = 
— 11.15 in both cases. The initial conditions for each numerical run 
are n h (0) = %>(0) = and n s (0) = e — > 0. All other parameters 
used are the same as in Fig. [3] 

of hard photons on the external ones must be taken into 
account. 

Thus, two more terms appear in the equations of hard 
and soft photons and the corresponding system now is (sys- 
tem S3): 

hp = Qpo — — — o- p7 n p n cx — <7 p7 n p n s (40) 

hh = — rih + An p n cx + An p n s — Chn s nh — C'^n^n^ (41) 
h s = — n s + C s n s n h + C' s n cx n h , (42) 
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Figure 9. Time evolution of the proton number density for pro- 
ton injection rate logQ po = —11.15 and two values of the ex- 
ternal number density, i.e. n cx = 2 (solid line) and n ex = 4 
(dashed line). The initial conditions for each numerical run are 
>%(0) = n p(") = an d n s(0) = e — > 0. All other parameters used 
are the same as in Fig. [3] 



where CL = " rr and CL = . These terms are derived 
using energy conservation considerations as in §3. The index 
'ex' is used to remind that the photon-photon cross section 
has a different value depending on whether the absorbing 
targets are the soft or the external photons. 

The additional coupling terms make now an analytical 
study as the one presented in §3 cumbersome. Thus, all the 
results presented in this section are derived after solving 
numerically the stiff set of equations S3. 

Our solutions indicate that these extra terms do not 
change, at least qualitatively, the validity of the results of the 
previous section. Increasing either Q po or n ex above a cer- 
tain value, forces the system to fall into a steady state rather 
than oscillate. As a first step, we compare two cases: (i) with 
and (ii) without the extra terms of absorption. In both cases 
we have used B = 0.7 G, 7 P = 2x 10 7 , logQ po = —11.15 and 
n cx = 2. What differs in the two cases are the external pho- 
ton energies, that are taken to be (i) e = lOOy^ 1 and (ii) 
e = 7 P ~ 1; ^7 L respectively. Although the system shows is os- 
cillatory in both cases, there is an important difference that 
is better displayed when the time evolution of the soft pho- 
ton distribution is plotted (see Fig. 1 1 0|) . At early times, the 
soft photon number density evolves in the same way in both 
cases. However, in the first case the number of soft photons 
in the source does not reach a deep minimum as in the sec- 
ond case. This can be attributed to the existence of the ad- 
ditional linear injection term of soft photons, i.e., -r-Can ex nh. 
Since the latter depends only on one time- varying parameter 
(rih), it adds soft photons into the system at a non negligi- 
ble rate, in contrast to the non-linear term +C a n a rih, that 
depends quadratically on the time-varying densities. 

This is also reflected on the shape of the limit cycles 
in the plane n p — n s of the phase space, which now appear 
tighter than before. This is shown in Fig. 1111 

As already discussed in section 2.2, the existence of an 
initial external photon distribution makes the role of quench- 
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Figure 10. Time evolution of the soft photon distribution n B 
in the case where hard photons are being absorbed on both the 
external and the soft photons (solid line) and in the case where 
they are being absorbed only on the latter (dashed line). Time 
is measured with respect to t,, that corresponds to the end of a 
transient phase that the system goes through, before it settles to 
its periodic state. For reasons of clarity, the transient phase is not 
shown. For the parameters used, see text. 
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Figure 11. Two dimensional plane logn p — logn s of the phase 
space in the case where hard photons are being absorbed on both 
the external and the soft photons (solid line) and in the case 
where they are being absorbed only on the latter (dashed line). 
Same parameters used as in Fig. 1101 

ing less clear. However, the way the equations of system S3 
are written, allows us to study separately the linear and non- 
linear absorption of hard photons by artificially deleting the 
non-linear terms of quenching in the equations of hard and 
soft photons, i.e., — Chn s rih and +C s n s rih respectively. For 
this purpose we examine two cases, that differ only at the 
proton injection rate. Figure [12] shows the hard photon 
compactness as a function of time for both cases, with panel 
(a) corresponding to the case with the larger proton injec- 
tion rate. The horizontal dotted line corresponds to , the 
solid line shows £h when both absorption channels operate, 
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whereas the dashed line shows £h when we ignore the soft 
photons produced due to non-linear quenching. Comparison 
of the solid and dashed lines leads to the conclusion that the 
automatic quenching of hard photons becomes dominant in 
the dynamics of the system from the instant that £h > £^ . 
Thus, even if quenching cannot be distinguished from the 
linear absorption when both operate, still remains an 
intrinsic property of the system. 

Another conclusion drawn from Fig. [T2] is that even in 
the absence of automatic quenching the system can exhibit 
damped oscillations. These example cases show clearly that 
the large-period limit cycle behaviour exhibited when au- 
tomatic quenching operates in parallel to the linear one, is 
replaced by an expontential growth that saturates, in the 
case where quenching is artificially omitted. On the other 
hand, the small-period limit cycle behaviour is replaced by 
damped oscillations of small amplitude. In other words, the 
combination of linear and non-linear absorption of 7-rays 
seems to intensify the temporal variability of the system. 
Finally, the absorption of hard photons is more efficient in 
the case where both the linear and non-linear channels of 
absorption operate than in the case where hard photons are 
being absorbed only on the external photon population. This 
can be deduced from the fact that, in the former case the 
averaged €h over a period is suppressed by at least one or- 
der of magnitude (see solid and dashed lines in panel (a) of 
Fig.II2). 

Two operating regimes of the system, each of them hav- 
ing its own properties, can be defined, depending on the 
relative strength of the two absorbing channels: (i) a linear 
and (ii) a non-linear one. While being in the linear regime, 
£h grows at first exponentially and eventually saturates. The 
transition from the linear to the non-linear operating regime 
depends on the external number density or equivalently on 
the optical depth r cx for absorption of hard photons on the 
external ones. In our analysis, t ox is simply given by 

r cx = o- 77 (e h eo)n ox (e ). (43) 

For small optical depths, i.e., r cx <C 1, this transition is 
abrupt in the sense that the system changes its tempo- 
ral behaviour completely. From the instant that automatic 
quenching becomes the dominant channel for absorption, the 
system exhibits limit cycles of large period and amplitude. 
On the other hand, for large optical depths, i.e., r cx > 1, the 
transition is smooth, since the system shows no limit-cycle 
behaviour. The non-linearity in the temporal behaviour be- 
comes evident by damped oscillations that reach a steady 
state in a few crossing times. Thus, the flaring behaviour of 
the system can be severely suppressed whenever the opti- 
cal depth for absorbing 7-rays on external photons is large. 
Panels (a) and (b) of Fig. [T3] show £^ as a function of time 
for two cases, where r cx = 1.3 and 0.13 respectively. In each 
panel, lines of different type denote different proton injection 
rates. The lightcurves depicted with dashed-dotted lines in 
both panels are obtained while the system operates in its 
linear regime. The dotted lightcurves exemplify the transi- 
tion to non-linearity, which is more abrupt for the case in 
panel (a) than the corresponding in panel (b). 

A new feature that appears through the study of sys- 
tem S3 is the dependence of the period T, if this exists, on 
the energy of the external photons e . This is exemplified in 
Fig. 1141 where the period of the oscillations varies with e 
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Figure 12. Time evolution of hard photon compactness when (i) 
7-rays are being absorbed both on external and on automatically 
produced soft photons (solid lines) and when (ii) the non-linear 
terms of absoprtion are artificially omitted (dashed lines). The 
dotted line in both panels corresponds to . The cases presented 
in panels (a) and (b) differ only at the proton injection rate, which 
is taken to be Q po = 3.2xl0 -11 and Q po = 4xl0 — 12 respectively. 
The rest of the parameters used in both panels are the same and 
equal to: B = 0.75 G, 7 P = 2.65 X 10 7 , e a = 1CT 5 and n cx = 1. 

almost like T oc l/er 77 (e ) - see eq. (|35|) . The minimum pe- 
riod is found at e that corresponds to the maximum value 
of the cross section for photon-photon absorption; for this 
value the absorption of hard photons becomes most effec- 
tive. The dashed line shows the dependence of the period 
on e for a higher density of external photons. In this case a 
gap appears for values of e that correspond to high values of 
the cross section around its peak. The evolution of the sys- 
tem there, is characterized by damped oscillations that lead 
eventually to a steady state. This result is to be expected 
within our analysis of §3. We remind that the system passes 
through well-defined stages, as one of the parameters Q po 
or rt cx increases: oscillations with large period — > oscillations 
with small period — > damped oscillations leading to a steady 
state. 

The fact that we find a clear analogy between the pe- 
riod and the inverse of the cross section for photon-photon 
absorption is a direct consequence of the simplifications we 
have made in the problem so far. However we note that, if 
were to relax our assumptions, i.e., use full expressions for 
the cross sections and emissivities, and treat the problem nu- 
merically, we would still have retained the basic conclusions 
of Fig. [HI 
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Figure 13. Hard photon compactness as a function of time for 
two cases with t cx = 1.3 (panel a) and 0.13 (panel b). In both 
panels the different light curves are obtained by increasing the 
proton injection rate. The transition from the linear to the non- 
linear operating regime of the system is clearly seen. Specifically, 
for panel (a) we have used: Q po = 10" 11 (dash-dotted line), 2 X 
10" 11 (dotted line), 4 X 1CT 11 (dashed line) and 8 X lCT 11 (solid 
line). The corresponding values for panel (b) are: Q po = 10~ 12 
(dash-dotted line), 1.6 X 10~ 12 (dotted line), 4 X 10~ 12 (dashed 
line) and 2.5 X 10 — 11 (solid lines). The rest of the parameters 
used are the same for both panels: B = 0.75 G, 7 P = 2.65 X 10' 
and to = 10~ 3 . 



We have also found that the period of the limit cycles 
depends not only on e but also on other parameters, that 
affect the value of the cross section of photon-photon absorp- 
tion even indirectly, as the magnetic field strength or/and 
the proton energy (see eqs. © and (©). The trend is 
the same as the one shown in Fig. 1141 where e in the hori- 
zontal axis should be replaced by the corresponding varying 
parameter. 



4.2 Inverse Compton scattering 

Thus far, we have assumed that the created pairs from 
photon-photon absorption act as an agent, transfering the 
energy from hard to soft photons through synchrotron radi- 
ation. However, if the compactness of soft photons becomes 
comparable or larger than the compactness of the magnetic 
field, i.e., l s > Ib, then there are two cooling channels for the 
secondary electrons: (i) the 'synchrotron' one that results to 
the production of soft photons e s and (ii) the 'inverse Comp- 
ton' (ICS) one that results to the production of high energy 




Figure 14. Dependence of the period T on the energy of the 
external photons e for number densities n cx = 1 (solid line) 
and n cx = 3 (dotted line). The rest of the parameters used are: 
Q po = 10" 10 , B = 3.57 G, 7 P = 9 X 10 6 and e Q = 8 X 10" 5 . 



photons ei cs - note that in general ei cs 7^ eh. Thus, the en- 
ergy lost by hard photons is only partially injected to the 
soft photon population. Because of this, the constant C s of 
the injection term in eq. (|42[) should be replaced by 



<?B + 34(l + 4e s7 e 



-3/2 ■ 



(44) 



where 7 C = eh/2 and the multiplication factor of l B ac- 
counts approxi mately for the Klein-N ishina cutoff effect up 
to e s 7 c < 10 4 (Modcrski ct al. 2005). Inspection of system 
S3 together with the expression (144 \ shows that the inclu- 
sion of ICS adds more non-linear terms to the problem. 

It is beyond the scope of the present work to proceed 
to a semi-analytical study of the above system. However, it 
is worth mentioning some qualitative effects of ICS on the 
dynamics of the system. In general, ICS acts as a damping 
term for soft photons whenever £ s > £b ■ 

Let us assume first, that we artificially switch-off ICS, 
and find a set of parameters that lead our system to a limit 
cycle behaviour as discussed in the previous section. If we 
keep the parameters fixed to these values and switch-on ICS, 
then there are three possible ways for the evolution of the 
system : 

(i) The limit cycle behaviour remains, although the sys- 
tem oscillates with a smaller period. 

(ii) The limit cycle behaviour is damped and the system 
finds its steady state after a number of oscillations. 

(iii) The system falls quickly into a steady state before 
showing any oscillations. 

The resulting behaviour of the system depends on the ra- 
tio £ s /£b and on whether or not the scatterings occur in 
the Klein-Nishina regime. Figure [15] exemplifies the above 
remarks. The solutions shown in Fig |15l are obtained after 
integrating the system of equations S3 and incorporating 
ICS in the approximate way described in this section. In 
both cases, the system starts with £ s <C £b but eventually 
reaches a state where £ s > Ib. The difference between the 
cases above is the parameter a;i cs = e s7 o, that denotes how 
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Figure 15. Soft photon density n s as a function of time for proton 
Lorentz factors 7 P = 2.1 X 10 7 and 2.65 X 10 7 (panels a and b 
respectively), derived for two cases: (i) ICS is artificially switched- 
off (solid lines) and (ii) ICS is taken approximately into account 
(dashed lines). Other parameters used for this plots are: B = 0.75 
G, logQpo = —11.15, n cx = 2 and e Q 



- 1 



deep into the Klein-Nishina regime the scatterings occur. 
Cases presented in panels (a) and (b) correspond to values 
a;ics = 7.6 and 87 respectively. In the former case, the damp- 
ing effect of ICS is evident, whereas in the latter case the 
system's evolution is not much affected because of the sup- 
pression of the scatterings. A small decrease in the period 
and in the amplitude of the oscillations is however evident. 
In the following section, where we treat numerically the full 
problem, we present a case that exemplifies the effects of 
ICS. 



5 NUMERICAL APPROACH 

All of our previous results were verified in an independent 
way by using the numerical code described in MK95 after 
selectively omitting various processes, as to make the code 
analogous to the systems described in the previous section. 

We proceed next to solve numerically the full system 
of eqs. (|15[) . Our aim is to present only some characteris- 
tic examples which will support our previous analysis and 
show also the effects of the processes we have neglected in 
our analytical approach, most notably of inverse Compton 
scattering, on the dynamical behaviour of the system. We 
leave a comprehensive search of the parameter space for a 
future paper, where the interplay between proton injection 
and external photons will be fully addressed. 

We have used the numerical code described in MK95 
and MPK05 which has been updated as to make use 
of the full rates for the secondary electron-positron pairs 
and photon production in photopion interac tions as ob- 
tained from the SOPHIA Monte Carlo code (|Miicke et al.1 



2000). Details on these will appear el sewhere - Dimitrak- 
oudis et al. in pr eparation; see also iDimitrakoudis et al.l 
|to appear in 2012 ). Therefore the updated version of the 
code can treat accurately the two major hadronic processes, 
i.e., photopair and photopion, in addition to the leptonic 
ones. Given the difficulty that these two processes pose in 
modelling, we consider this as a major improvement. 

We solve therefore three coupled equations, for protons, 
electrons and photons including all relevant processes be- 
tween the three species - note that in the numerical code 
there is no need to treat the hard and soft photons through 
separate equations. Another difference with the analytical 
treatment is that in place of the external photons we use the 
photons produced by the proton synchrotron radiation. This 
was done because this process can produce the seed photons 
self-consistently without the need of introducing more free 
parameters. Furthermore, the proton synchrotron radiation, 
for magnetic fields of ~ 10 G and protons with Lorentz fac- 
tors 7 P > 10 6 is produced mainly in the soft energy range of 
the spectrum and cannot/should not be neglected. 

Fig. [16] shows four cases that differ only in the pro- 
ton injection rate. Thus for panels (b), (c) and (d) Q po was 
increased by a factor of 2, 3 and 4 respectively over its cor- 
responding value of panel (a). The latter one was chosen in 
such a way as to make the system exhibit large period limit 
cycles. The period starts decreasing with increasing Q po - 
as a matter of fact the period gets exactly to a half of its 
previous value as Q po is increased by a factor of 2. This 
behaviour degenerates into a damped oscillation - steady 
state mode with increasing Q po (panels c and d). This is 
exactly the behaviour we found in our analytical treatment 
- see Fig. [5] Therefore, despite the plethora of the physical 
processes introduced, the feedback loop still operates. 

For low values of the magnetic field, ICS acts as a fric- 
tion mechanism stabilizing the system and letting it reach 
quickly steady state. Fig. ll7l shows an example where we ran 
the code for two cases. In the first, all processes were taken 
into account (dashed line), whereas in the second case, ICS 
was artificially switched-off (solid line). In the former case 
the system, after an initial peak, falls quickly into a steady 
state. In the latter it shows the limit cycle behaviour. The 
two cases are identical until the time when the compact- 
ness of soft photons becomes large and cannot be neglected 
any further in the electron cooling. This occurs around the 
time that the first peak in appears, where the soft photon 
compactness becomes by a factor of ~ 6 x 10 3 higher than 
the magnetic one. It is worth comparing the result shown in 
Fig. [T7] with the one plotted in panel (a) of Fig. 1151 where 
ICS was taken into account in an approximate way. In both 
cases the qualitative results are the same. If we use a higher 
value of the magnetic field for the example case shown in 
Fig. 1171 we find that ICS alters the periodic behaviour only 
by decreasing the period. The temporal behaviour we find 
in this case, can be very well described by the correspond- 
ing one shown in panel (b) of Fig. 1151 Suppresion of ICS 
because of Klein-Nishina cutoff effects does not play an im- 
portant role in this case, since we do not assume monoener- 
getic electron and photon distributions. Thus, most of the 
scatterings occur in the Thomson regime and Klein-Nishina 
cutoff effects are now small corrections. Finally, we note that 
the limit cycle behaviour of the system shown in Fig. [16] re- 
mains intact, although ICS was taken into account, since for 
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Figure 16. Time evolution of the system for four different pro- 
ton injection rates or equivalently proton compactnesses, starting 
with 4 nj = 4.7 X 10~ 5 in panel (a). In panels (b) to (d) £ p n > is in- 
creased over its previous value by an integer multiple of its initial 
value. The rest of the parameters used for this plot are: B = 10 G, 
R = 10 16 cm and 7 P = 3 X 10 7 . Solid and dashed lines show d p 
and £j respectively, where d p = f d7p7 P np(7 P ). 



the parameters used, the magnetic energy density is always 
larger than the soft photon one. 



6 RELEVANCE TO ASTRO-PHYSICAL 
SOURCES 

We turn next to examine the ideas presented in the previous 
sections in the context of possible applications of astrophysi- 
cal interest. Whenever the system operates in the subcritical 
regime or the 7-rays are being absorbed mainly through the 
linear absorbing channel (see §4.1), our steady state spectra 
are similar to those presented in the literature - for example 
sec iBottcherl l|2007l ). On the other hand, if the system be- 
comes supercritical and the absoprtion of 7-rays is strongly 
non-linear, then new possibilities open for astrophysical ap- 
plications. These, according to §3 and §4, can be summarized 
in the following: 

(i) All the components of the hadronic system exhibit an 
intrinsic variability with a well denned period, although the 
source is stationary. 

(ii) The system reaches a steady state after going through 
a damping oscillatory phase. In this case a soft photon com- 
ponent emerges since a significant fraction of the energy 
stored in protons is transferred to lower energy photons via 
quenching of the 7-rays. At the same time the hard photon 
compactness reaches a limiting maximum value. 



0- 
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Figure 17. Photon compactness £j as a function of time for two 
cases where, (i) inverse Compton scattering is taken into account 
(solid line) and (ii) it is artificially switched-off. The parameters 
used for this plot are: B = 3.2 G, R = 10 16 cm and 7p = 3 X 10 7 , 

e p ni = 7 x icr 5 . 



Both cases can, in principle, have relevance to astro- 
physical sources emitting in high energies, such as AGN. 
Time variability is a denning property of blazars and, 
in most cases, i t shows a quite comp l ex pa ttern (e.g., 
iMukheriee et all Il999l : lAharonian et al.l l2007h . We note 
that, even if the observations seem to contradict our results 
(see point (i) above), we have found that even small ampli- 
tude variations of the proton injection rate can lead to much 
more complicated time profiles than the ones presented so 
far. This is a promising topic. Therefore, a detailed study of 
the system towards this direction is required and is going to 
be the subject of a future work. 

In the rest of this section we will focus on the second 
point. We will show specifically how we can apply our results 
in order to set an upper limit to parameter values used in 
the modelling of AGN multiwavelength spectra. This can be 
seen as an extension of the application presented in PM11. 

The luminous blazar 3C 279 is a good example. A re- 
cent comprehensive review of observa tions can be found in 
iBottcher. Reimer. fc Marscherl |2009l '). Here we will mainly 
focus on the 2006 campaign, which discovered the source at 
VHE 7-rays, showing a high TeV flux (Albert et al. 2008), 
while the X-rays were at a much lower level. 

Let us consider a spherical source of radius R moving 
with a Doppler factor 5 with respect to us and containing a 
magnetic field of strength B. We further assume that ultra- 
relativistic protons with a power law distribution of index 
s are being constantly injected into the source with a rate 
given by 

Qp = Qpo7p S H{~f p - 7mm)-ff(7ma X - 7p), (45) 

where 7 m i n and 7 max are the lower and upper limits of the 
injected distribution respectively. Q po is the normalization 
constant and is also directly related to the proton injection 
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Figure 18. Multiwavelength spectra of 3C 279 obtained in the 
context of a pure hadronic model for R = 3 X 10 16 cm, B = 40 
G, S = 20, 7 min = 5 X 10 s , 7 m ax = 5 X 10 9 , s = 2.2 and three 
proton injection compactnesses: t p s = 10~ 5 (solid line), 2 X 10~ 5 
(dashed line) and 4 X 10 -5 (dotted line). The symbols represent 
the observational data from February 2006. 

compactness as: 

if = Q po m p c 2 ^ 7 ^ + 2 2 J g 7 "- (46) 
or in terms of dimensionless quantities 

t n 3 = l Qpo 7 m ,n^ 7max _ (47) 

The transition of the system from the sub - to the super- 
critical regime can be better seen if initially there are no soft 
photons present in the source. For this, we try to fit only the 
TeV emission by considering a narrow power law proton en- 
ergy distribution and no primary electron population. Thus, 
the multiwavelength spectrum is purely the result of proton 
primary and secondary emission. 

Figure [18] shows the multiwavelength spectra obtained 
using the numerical code described in §5 for R = 3x 10 16 cm, 
B = 40 G, 5 = 20, 7 min = 5 x 10 8 , 7max = 5 x 10 9 , s = 2.2 
and for three values of the proton injection compactness 
starting with t p s = 10 -5 (solid line) and increasing it over 
each previous value by a factor of two. For this set of param- 
eters the produced 7-rays lie in the GeV-TeV regime. The 
spectrum for the lowest value of t p nl (solid line) is obtained 
while the system is subcritical, and it is the only one that 
can give an acceptable fit of the TeV emission while at the 
same time does not violate the optical and X-ray observa- 
tions. We see that an increase of t p 3 by just a factor of two 
drives the system into the supercritical regime. The onset 
of supercriticality is acompanied by the emergence of a soft 
emission component (dashed line), that becomes dominant 
for an even higher t p 3 (dotted line). The over-production of 
soft photons even in the second case is excluded directly by 
the observations, setting an upper limit to the proton injec- 
tion compactness (for the specific example, ^p" J max = 10 -5 ). 

If we were to use a lower value of 5 to obtain the fit 
in the TeV energy range while keeping fixed the magnetic 
field strength, we would require a higher value of £ p nl , since 
L bs oc (5 4 Li nt , where L b s and Li nt are the luminosities in 



the observer's and in the comoving frame respectively. This 
choise of parameters would drive the system deep into the 
supercritical regime violating the optical and X-ray observa- 
tions. PM11 have used similar arguments to set constraints 
on the Doppler factor <5 by using an ad-hoc 7-ray injected 
luminosity. Here we go one step further since the injected 7- 
rays are related to a physical production mechanism. Thus, 
in this case we can set a limit on both 5 and £ p nl . Moreover a 
potential flaring event observed only in the GeV part of the 
spectrum could not be fitted by just increasing the proton 
injection luminosity, since this increase would affect also the 
optical and X-ray part of the spectrum, as the example of 
Fig. [T8l suggests. 

Thus, the effects of the underlying feedback mechanism 
can prove to be useful in setting lower limits to parameters, 
such as the Doppler factor 8. A systematic search of the 
parameter space is however out of the scope of the present 
paper. 



7 DISCUSSION 

Hadronic models have been extensively used for explaining 
AGN non-thermal emission while recently thay have been 
applied, as well, to other compact objects. An interesting, 
but overlooked property of hadronic systems is their dynam- 
ical behaviour, which results from some underlying feedback 
mechanisms. In the present paper we have isolated and stud- 
ied analytically one such loop that involves proton-photon 
pion production and photon quenching (SK07; PMU) of 
the produced 7-rays. Gamma-ray quenching results in au- 
tomatic production of soft photons which then feed back 
producing more proton-photon cooling via pion production. 

We have remarked that if protons are considered sta- 
tionary in the source - an often made assumption, there 
are parameter space regimes which are characterized by an 
exponential growth of internally produced photons, making 
the system inherently unstable since the condition for proton 
stationarity is violated by the losses caused by the runaway 
photons - here the analogies to the 'Compton catastrophe' of 
leptonic plasmas are evident. The kinetic equation approach, 
which allows proton cooling to be explicitly taken into ac- 
count, is suitable for studying in this case, the properties of 
the system. Protons, secondary electrons and photons, i.e., 
the main three components of the system, can be described 
by three coupled partial integro-differential equations. This 
'kinetic equation treatment' has many advantages, as it is 
both energy conserving and time dependent. 

In order to simplify the system of equations and make 
an analytical treatment possible we made a number of as- 
sumptions. As a first step we have retained only simplified 
expressions of the key processes by using <5-functions for 
the different particle distributions appearing in the prob- 
lem. Furthermore, we used approximate expressions for the 
cross sections - see eqs. (0 and (|35|) . However, one of the 
major simplifying assumptions we made was the elimination 
of the electron kinetic equation from the system of eqs. (| 15f) . 
The rationale for this is that the electron cooling timescale, 
for typical values of the magnetic field and electron Lorentz 
factors, is much smaller than the crossing time of the source. 
Thus, electron cooling is considered to be instantaneous. 

Ignoring absorption of the hard photons on the external 
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ones, we have found that an increase of the proton injection 
rate leads to an analogous increase of the proton density and 
of the hard photon luminosity resulting from photopion. If 
the latter does not get to a high value as to trigger quench- 
ing, the system is linear and reaches a steady state. If, how- 
ever, the combination of the initial parameters is such that 
leads the hard photons past the quenching threshold, then 
the latter are automatically absorbed and the soft photons 
which are spontaneously produced serve as targets for extra 
proton cooling. In this case we have shown analytically by 
performing an eigenvector/eigenvalue analysis, that for hard 
photon densities close but above critical the system goes 
through a limit cycle behaviour of the prey-predator type. 
For even higher hard photon densities, the system reaches 
a steady state which is achieved after protons and photons 
exhibit a series of damped oscillations. This steady state oc- 
curs at very different values from the ones achieved while the 
system is subcritical and this discontinuity is another indica- 
tion of the system's supercriticality. It is interesting to note 
th at duty-cycle behaviour in hadronic systems was reported 
bv lStern fc Svensso 3 <ll99lT> using a Monte Carlo code and 
by MPK05 using a kinetic equation approach. However both 
papers were numerical and the authors, while giving ample 
physical reasoning for the behaviour, fell short of presenting 
a mathematical proof. Here for the first time we present such 
an interpretation and show beyond doubt, that hadronic sys- 
tems can indeed exhibit the aforementioned behaviour. 

As a next step we have shown semi-analytically that 
when absorption of the produced hard photons on the ex- 
ternal ones is included, then the behaviour of the system de- 
pends also on the corresponding optical depth. If this takes 
low values, then the system behaves as described above since 
the non-linear processes continue to play a dominant role in 
the dynamics of the system. If, on the other hand, the opti- 
cal depth takes high values, the systematic depletion of hard 
photons tends to stabilize the system which reaches a steady 
state after it goes through a damped oscillation mode, i.e., 
no limit cycle behaviour was found in this case. It is in- 
teresting however to note that, even in this case, the hard 
photons cannot reach a steady state above the critical value 
of quenching. This can only mean that quenching remains a 
fundamental intrinsic property of the system. 

Furthermore, if we are to add more physical processes 
to the system, these tend to stabilize it as they redistribute 
part of the radiated energy away from the operating feed- 
back loops. For example, inverse Compton scattering can 
act as a competing energy loss mechanism for electrons 
to synchrotron radiation. If it dominates, then the system 
changes behaviour moving faster to a steady state. However 
for strong enough B-fields, we found that the system retains 
the analytically derived properties. 

The above results, were also verified by using a numer- 
ical code where the full expressions for the emissivities of 
the various radiative processes and for the various relevant 
cross sections were used. While the details change from our 
simplified analytical approach, we were able to confirm qual- 
itatively our analytical results which predict the transition 
from the subcritical linear regime to the supercritical oscil- 
latory one with increasing proton density. In addition, we 
were able to verify the role of other processes, like inverse 
Compton scattering, which we had taken in our analytical 
treatment only approximately into account. We note that 



the qualitative analogies between the results of the two treat- 
ments justify also a posteriori and in an independent way 
the validity of our simplifying assumptions. 

Our results indicate that higher proton injection rates 
tend to push the system into the non-linear regime with 
the external photons acting more as catalysts; on the other 
hand, higher external photon densities act on the opposite 
direction and tend to linearize the system. Preliminary nu- 
merical calculations (Dimitrakoudis et al. - in preparation) 
show that this trend remains, if one is to replace the <5- 
function distributions used in the present treatment with 
the more astrophysically relevant power-laws. 

Finally, as an example of astrophysical interest, we have 
used the numerical code described in §5, which can treat 
self-consistently both the non-linear development of EM cas- 
cades and proton cooling, in order to make a fit to the TeV 
emission of blazar 3C 279. We have shown that acceptable 
fits can be obtained only for high values of the Doppler fac- 
to 1 ' (<5 ft 20), given that typical values for the magnetic field 
in the context of hadronic models are considered (B ~ 10 
- 60 G). However, this is only an indicative example of the 
potential applications of the model. Obviously one needs to 
thoroughly study the parameter space before providing ex- 
act fitting values. Another potential direction is the study 
of the inherent variability signatures of the system in cases 
where the source itself is variable. 

The supercriticality related to the feedback mechanism 
studied in the present paper is by no means the only one that 
can occur in hadronic systems. If one was to replace the pho- 
topion interactions by another production mechanism of j- 
rays, e.g., proton synchrotron radiation, and study the same 
feedback loop outlined in §1, one would again find that the 
system enters a supercritical regime. However, the parame- 
ter values that would enable the transition from the subcrit- 
ical to the supercritical regime, as well as the the details con- 
cerning the transition itself, would differ from those shown in 
the present work, due to diferrent cross sections, emissivities 
and energy threschold criteria. Furthermore, we cannot ex- 
clude that even other loops operate as well in a hadronic sys- 
tem - see iKirk fc Mastichiadis! (1 19921 ). iDimitrakoudis et al.l 
(|to appear in 20121 ). Actually, more than one of the different 
feedback mechanisms can operate simultaneously in a 'real' 
hadronic system. Because of this, a comprehensive search of 
its parameter space demands the use of the numerical code 
described in §5 and it is going to be the subject of a future 
work. At any rate, we have shown that hadronic models con- 
stitute one more example in the growing list of dynamical 
systems and as such they need to be further investigated. 
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APPENDIX A: STABILITY ANALYSIS OF THE 
TRIVIAL STATIONARY SOLUTION 



Consider the system (S) 

—x + An cx z + Azy — Chl/x 
-y + C s yx 



x = 

y = 

i = 
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zy 
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where variables (x,y,z) stand for (nh,n s ,n p ). The fixed 
points of the system, P%{xq ,y$ ,z$ ) with i — 1,2,3 can 
be found by setting x = y — z — 0. The first fixed point 
is nothing more than the trivial steady state solution of the 
system in the case of no quenching: 
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We make the convention that y is the positive real root 
that has physical meaning. It is interesting to examine the 
behaviour of the system when it is slightly perturbed by its 
steady state, i.e., x = x^ + x',y = y',z — Zq + z' . For 
this we linearize system (S) with respect to the perturbed 
quantities: 

x> = -x' + (Az\? - + An^z (A8) 

y> = (-l + Csa^V (A9) 
z' = -ff^z^V - G p z. (A10) 

For the second equation that is not coupled to the other 
two we find an exponential growth or decay y'(r) = j/'(0)e ST 
depending on the sign of s = — 1 + C s Xq . First suppose 
that s < 0. y — > holds for r > 1/s and the three lin- 
earized equations degenerate to two. The matrix of the two 
dimensional system is 



Mi = 







An cy 

— Gr> 



APPENDIX B: STABILITY ANALYSIS OF THE 
NON-TRIVIAL STATIONARY SOLUTION 

In general, the topology of a vector field near its fixed points 
can be studied through the Jacobian matrix of the vec- 
tor field at the corresponding points. Specifically, the clas- 
sification of fixed points in different types is made by an 
eigenvalue/eigenvector analysis of the Jacobian matrix. This 
analysis is widely used for two-dimensional vector fields, 
leading to a few types of fixed points. This is not the case 
for three dimensional systems, where the classification of the 
fixed points in types is more complicated. In our work we 
have adopted the classification presented in iTheisel et al.l 
{200S) . 

Here we apply the eigenvalue analysis to the three di- 
mensional vector field v = (nh,n s ,w p ) T . We consider the 
second non-trivial stationary solution of the system S2. 
Thus, the set of equations (S2) after linearization at the 
point P2 can be written in the form 
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= M 2 


: 









where the matrix M2 is given by 
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The eigenvalues are the roots of its characteristic polynomial 

(Bl) 



P(X) = A 3 + aiA 2 + a 2 A + a 3 



where 
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The number of real and complex roots of the equation 
P(A) = can be determined, if one knows the signs of the 
constants. For values relevant to our physical problem, one 
finds that (i) either all three constants are positive or (ii) 
only a-i is negative. In the first case the polynomial has 3 
negative real roots or 1 negative real root and two complex 
conjugates, while in the second case there is always 1 neg- 
ative real root and 2 complex conjugates with positive real 
parts. 

Figure [BTJ shows the dependence of constants a, on Q po 
and n ex . Since a logarithmic scale is used, the negative values 
of a,2 are not shown. 



with determinant A(Mi) = G p > and trace Tr(M 1 ) = 
— 1 — G p < 0. Thus, in this case the point Pi is stable. This 
simple analysis does not apply in the case of s > 0. However, 
it can be easily shown that both z' and x' are oc e ST for large 
enough times. Thus, in this case all the perturbed quantities 
grow with time and Pi can be characterized as unstable. 



APPENDIX C: DERIVATION OF THE 
CRITICAL 7- RAY COMPACTNESS 

Let us assume that hard photons are being injected into a 
spherical source with a constant rate Q™ 1 that corresponds 
to a compactness t^ s and that no soft photons are initially 
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Figure Bl. Dependence of the constants of the characteristic 
polynomial P(X) on Q po for n ex = 2 (top panel) and on n e x for 
logQpo = —11.15 (bottom panel), for a range of values relevant 
to the physical problem. In both panels solid, dashed and dotted 
lines represent the constants cii, (12 and 123 respectively. 



present in the source. Automatic quenching of hard pho- 
tons is possible if the injected compactness exceeds a certain 
value. In this case, electron-positron pairs are being created 
spontaneously in the source, emitting synchrotron radiation. 
Hard photons then undergo further absorption on the afore- 
mentioned soft photons and a non-linear loop of processes 
begins to operate. The equations that describe the above 
physical system can be written in the following form: 



rih = Qh° J - rih - C h n s n h 
n s = —n s + C a n a nh, 



(CI) 
(C2) 



where the constants Ch and C s are denned in eq. (|33l) . 

There is a trivial stationary solution of the system (|C2|) : 
fih = Q'h 3 , n s = 0, that corresponds to the free propaga- 
tion of hard photons through the source, where pairs and 
soft photons are absent. The Jacobian matrix evaluated for 
this solution has two real eigenvalues. For Q™ J < 1/C S both 
are negative. Thus, the solution is stable. Moreover, in this 
case the system has no other physically acceptable solution, 
i.e., n s > 0. However, if QJJ 1J > 1/C S one of the eigenval- 
ues becomes positive and the solution with a zero soft pho- 
ton population becomes unstable. Even a perturbation in 
the initially absent soft photon distribution is sufficient for 



its subsequent growth. In this region, a second stationary 
solution of the system (|C2[) appears. This is fih = 1/C S , 
n B — (q'^' — \ /Chfih and one can show that for this re- 
gion both eigenvalues of the corresponding Jacobian matrix 
are negative, i.e., the solution is stable. 

Summarizing, the critical injection rate is 1/C S , which 
can be transformed to the critical compactness: 



3a 



(C3) 
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